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Abstract 

We study reduced-order models of three-dimensional perturbations in linearized channel flow using bal- 
anced proper orthogonal decomposition (BPOD). The models are obtained from three-dimensional sim- 
ulations in physical space as opposed to the traditional single-wavenumber approach, and are therefore 
better able to capture the effects of localized disturbances or localized actuators. In order to assess the 
performance of the models, we consider the impulse response and frequency response, and variation of the 
Reynolds number as a model parameter We show that the BPOD procedure yields models that capture the 
transient growth well at a low order, whereas standard POD does not capture the growth unless a consider- 
ably larger number of modes is included, and even then can be inaccurate. In the case of a localized actuator, 
we show that POD modes which are not energetically significant can be very important for capturing the 
energy growth. In addition, a comparison of the subspaces resulting from the two methods suggests that 
the use of a non-orthogonal projection with adjoint modes is most likely the main reason for the superior 
performance of BPOD. We also demonstrate that for single-wavenumber perturbations, low-order BPOD 
models reproduce the dominant eigenvalues of the full system better than POD models of the same order. 
These features indicate that the simple, yet accurate BPOD models are a good candidate for developing 
model-based controllers for channel flow. 
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I. INTRODUCTION 



Many techniques for developing practical controllers for fluids require models of the system 
that are both tractable and that describe accurately the flow physics for the given flow regime. 
One of the problems of great interest in flow control is drag reduction in shear flows. Drag in- 
creases drastically as flow transitions from laminar to turbulent, making turbulence suppression 
or inhibition of transition to turbulence promising strategies for drag reduction. Both open-loop 
and closed-loop strategies have been used in control of channel flow. Recently, Min et al ll22ll 
have achieved sub-laminar drag through open-loop control using a traveling wave actuation on the 
channel walls. Lee et al [|20ll have reported drag reduction using closed-loop controllers developed 
for the linearized flow, while Joshi et al ifTSl and Hogberg et al [fT2| have demonstrated significant 
reduction in the energy of the perturbations to laminar flow using closed-loop controllers. Despite 
these recent successes, the adequate modeling of the fluid flow and the actuators in a framework 
useful for practical control design is still a challenge. In order to use model-based control strate- 
gies, one needs an accurate description of the system dynamics and the actuation, from a model 
of sufficiently low order to allow practical implementation. The goal of this paper is to explore 
improved techniques for developing such reduced-order models, in the context of a transitional 
channel flow. 

The POD/Galerkin method has been used extensively for reduced-order modeling of fluid prob- 
lems lfT3ll23ll3n[34l[35ll . In this method, an empirical basis of orthonormal eigenfunctions is ob- 
tained from experimental or simulation data, and the Navier-Stokes equations are projected onto 
this basis. For fluid problems, this basis is optimal in terms of capturing the energy of the flow: the 
most significant modes are the ones that carry most of the kinetic energy. Although this method 
is applicable to many different types of flows, and computationally tractable for very large data 
sets, it can result in inaccurate low-order models since the most energetic modes are not always 
the most dynamically significant ones, as for example in the case of acoustic modes in cavity 
oscillations Bill . Improved POD models can in some cases be obtained by a careful selection of 
modes [[35l . removal of symmetry from the data using Fourier or traveling modes [[T4l[32ll or inclu- 
sion of shift modes [|25l . but these techniques are often ad hoc and typically require extensive fine 
tuning. In addition, there is no standard method for incorporating actuation into POD models, and 
if using these models for control design, one needs to be careful to respect the region of validity 
of the model in the presence of actuation [36|. 
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On the other hand, many tools from the control theory community have been developed for 
model reduction of linearized dynamics, often with a priori error bounds [|26l . Here we focus on 
one of these methods, balanced truncation, introduced by Moore [|24l|. and described in detail in 
standard references ([|4l|26l[39l). Balanced truncation involves transforming a state-space system 
to a coordinate system where the states that respond most strongly to inputs (most controllable 
states) are also the states that have the most influence on future outputs (most observable states). 
In this balanced realization, weakly observable and controllable states can be truncated to form 
reduced-order models that capture well the input-output behavior of the system. The required 
coordinate transformation is known as the balancing transformation. 

The main difficulty with using balanced truncation for large systems is computational expense. 
In a typical fluids problem, the number of states is on the order of 10^ or more, and since com- 
puting balancing transformations involves solving Lyapunov equations for matrices of dimension 
n X n (where n > 10^ is the number of states) as well as correspondingly large eigenvalue prob- 
lems, traditional approaches to balanced truncations may not be feasible. Recently, snapshot-based 
methods for computing balanced truncation have been suggested, both to extend the concepts to 
nonlinear systems [fTSl and to large linear systems [[30l . For systems with large numbers of out- 
puts (e.g., if the output is the entire state), the snapshot-based method requires a large number 
of adjoint simulations, which still makes the problem computationally intractable. A procedure 
we call Balanced POD (BPOD) [30], uses an output projection to reduce the number of neces- 
sary adjoint simulations, and has been shown to be a computationally feasible approximation to 
balanced truncation for examples of two- and three-dimensional perturbations to laminar channel 
flow 

Many recent works on transition in shear flows have focused on the large non-normal transient 
growth of exponentially stable linear perturbations to the laminar flow, which is thought to lead to 
the so-called 'subcritical' or 'bypass' transition [|2l|3l[6l[T6l|29l|33l|37l. A comprehensive treat- 
ment of the subject is given by Schmid and Henningson [|33l . Suppression of this large growth is 
of interest in control applications for drag reduction. Control and estimation of linearized channel 
flow was studied by Hogberg et al [|l2l. More recently, Akervik et al [fill studied control of a cavity 
flow using global eigenmodes. In the works by Farrell and loannou [7 J and Lee et al [|20l . balanced 
truncation was applied to linearized channel flow at particular wavenumbers where the standard 
algorithms are applicable, since the full system is one-dimensional. 

The contribution of this work is the application of BPOD to transitional channel flow for 
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localized disturbances, without modeling on a wavenumber-by-wavenumber basis. Although a 
lot has been learned from studying 'canonical' classes of perturbations (streamwise vortices, 
oblique waves, ToUmien-Schlichting waves, simple localized disturbances [[TOl ), the standard 
wavenumber-by-wavenumber analysis has its limitations when efficient implementation of closed- 
loop control is desired, since model reduction and control design would need to be performed at 
each wavenumber. The perturbations that arise in real flows often have complex three-dimensional 
structure with contributions at a wide range of wavenumbers as well as non-periodic or stochas- 
tic components. More importantly, for modeling of realistic actuation devices the wavenumber 
approach would be very complex except for the special case of actuators that excite only spe- 
cific wavenumbers. It is therefore of interest to have a method where disturbances with complex 
structure are modeled without Fourier decomposition, in particular for controls applications where 
physically reahzable localized actuators need to be considered. To the best of our knowledge, such 
low-order balanced truncation models have not been reported in literature. Although some suc- 
cess has been achieved for channel flow using control at all wavenumbers, the method we propose 
is able to extract dominant dynamics of the flow in physical space, resulting in simpler models. 
Moreover, this method for balanced truncation would be very advantageous for more complex 
geometries, such as spatially developing boundary layers. 

We first obtain models of a single- wavenumber perturbation computed by Farrell [3] in order 
to validate our numerical methods, and then apply BPOD to a localized body-force actuator. The 
most desirable features of a reduced-order model are close approximation of the dynamics of the 
original system, inclusion of actuation and validity over a wide range of parameters, and we show 
that our BPOD models have those features. We also show that BPOD models capture the effects 
of actuation better than standard POD and that the considerable improvement in the capturing of 
the dynamics for BPOD is due to the non-orthogonal projection used. 

The rest of this paper is organized as follows. In Section [III we briefly overview the two model 



reduction methods we use. In Section III we describe the governing equations and our choice of 



inner product for the adjoint simulations of the system. In Sections |IV] and |V} we present low-order 
models for a single-wavenumber optimal perturbation and a localized actuator and discuss their 
performance. Finally, in Section IVll we describe our conclusions and directions for future work. 
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II. MODEL REDUCTION VIA POD AND BALANCED POD 



A. Proper Orthogonal Decomposition (POD) 

We first give a brief overview of the POD/Galerkin method; details can be found in standard 
references lfT3l[34ll . The idea of Galerkin projection is, given a system 

X = fix), x{t) e X, (1) 

where X is a high-dimensional Hilbert space, to project onto a low-dimensional subspace S C X. 
Proper Orthogonal Decomposition determines an orthogonal basis for such a subspace, which is 
obtained by solving the eigenvalue problem for XX^ where X is a matrix whose columns are 
simulation snapshots x{tk) at some times tk. In this basis, we can represent the dynamics of x{t) 
as 

m 

x{t) = J2(^jitWj, (2) 

i=i 

where 9j are the time-independent basis functions (POD modes) and aj (t) are the corresponding 
time coefficients, which are obtained from 

d, = (%,/(r)). (3) 

A reduced-order model of order r can be obtained as a set of ODEs for the time evolution of these 
coeffecients by projecting the original system onto the most significant POD modes in the system 
(i.e. including only the first r modes where r < m). The method is applicable to both linear and 
nonlinear systems. For linear systems, the POD modes of data arising from the input-state impulse 
response are the most controllable modes of the linear system [|30ll . However, both controllability 
and observability are important for the input-output behavior of a system, and POD often fails to 
capture highly observable modes. On the other hand, balanced truncation does take into account 
both of these properties, and we describe this method next. 



B. Balanced POD 

Balanced truncation is a standard model reduction method flU |24l |39l used for stable linear 
input-output systems of the form 

X = Ax + Bu 

(4) 

y = Cx 
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where m G W = is the vector of inputs, y E y = W is the output, x E X = is the state 
vector (although in general all three spaces can be complex as well), and A, B, and C are matrices 
of appropriate dimension. The idea of balancing is to find a change of coordinates in which the 
controllability and observability Gramians, defined by 

/•oo /»oo 

W,= / e'^'BB-^e^^'dt, Wo = / e^'^'C+Ce^' dt, (5) 
Jo Jo 

are equal and diagonal. Here A~^, B^ and define the corresponding adjoint system. It should 

be noted that in general 7^ , the two being equal only when the inner product used to derive 

the adjoint does not have an associated weight. It can be shown that balanced truncation does not 

depend on the choice of the inner product on the state space X (see Appendix [A]), although it does 

depend on the choices of inner products for U and y . One then truncates the least controllable 

and observable modes, corresponding to the smallest eigenvalues of these Gramians. A detailed 

description of the Balanced POD procedure, which is a computationally tractable procedure for 

finding such a transformation, is given in [30J. In this method, one begins by computing snapshots 

of the impulse-state response of the system in Eq. (|4]) and the adjoint system 

z = A-^z + C+w (6) 

and stacking the direct and adjoint snapshots as columns of matrices X and Y (with appropriate 
quadrature weights [l30l '). One can show that the Gramians in Eq. (|5]) may then be approximated 
by empirical Gramians [fTSlI VFc,e and Wo^e, as 

^ iy,,e = Wo ~ Wo,e = YY+. (7) 

The key idea in the method of snapshots is to compute the transformation that balances the empir- 
ical Gramians (or at least the dominant directions of this transformation) without actually comput- 
ing the Gramians themselves, whose dimension is large. In this respect, the method of snapshots 
for BPOD [|30ll resembles the method of snapshots introduced by Sirovich ll34l for more efficient 
computation of POD modes. To compute the balancing transformation, one computes the singular 
value decomposition (SVD) of the matrix Y^X (see Appendix [a| for a discussion of Y^): 

Y+X = U^V^, (8) 

from which the balancing transformation $ and its inverse ^ are found by 
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The columns of $ are the balancing modes and the columns of \& are the adjoint modes, and the 
two sets of modes are biorthogonal. The entries of the diagonal matrix S are known as the Hankel 
singular values (HSVs). The non-orthogonal projection onto the basis of balancing modes using 
the adjoint modes is also known as Petrov-Galerkin projection. 

Note that a different procedure for approximating balancing transformations has also been used 
in ll38l . in which the Gramians are separately reduced (that is, low-rank approximations of Wc,e 
and Wo^e are first constructed, and then the balancing transformation for the rank-reduced Grami- 
ans is computed by an unspecified algorithm). However, this procedure is more computationally 
intensive than our procedure, and also gives worse results [301, since almost-uncontrollable modes 
may be strongly observable, so should not be truncated. 

C. Output projection 

If the number of outputs of the system is large, as in a typical fluids problem (where n = q, 
i.e. the output is the full state), the computation of the adjoint simulations of the system given by 
Eq. ([6]) may not be tractable, since one simulation is needed for each component of the output. A 
way to reduce the number of system outputs is to first project the output onto a low-dimensional 
subspace, i.e., taking y = PgCx, where Ps is an orthogonal projection onto a s-dimensional 
subspace of y, as suggested in [|30l . The system is now of the form: 

X = Ax + Bu 

(10) 

y = PsCx 

where s is the rank of the output projection. The projection Ps that minimizes the 2-norm of 
the difference between the original transfer function and the output-projected transfer function is 
given simply by the POD of the set of impulse-state responses [|30ll . This projection can be written 
as Ps = ©s©^, where columns of ©^ : M'* ^ 3^ are POD modes. Another way to write the system 
is as follows: 

X = Ax + Bu 

(11) 

y = etcx 

Here, the outputs of the system are just the coefficients of the POD modes of the system impulse 
response and y eW. This s-dimensional output carries the same information as the n-dimensional 
output y, which can be shown by Parseval's theorem. The corresponding adjoint system can now 
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be written as 

z = A+z + {QtC)+v 

(12) 

w = B^z. 

Note that if the output is the full state, so that C = I, and the adjoint is defined with respect to 
the standard L2 inner product, the initial conditions of the adjoint simulations are just the POD 
modes (columns of 6^). In practical computations, depending on the choice of inner product 
used in defining the adjoint system, and on the numerical quadrature method (for example, if the 
computations are done using Chebyshev polynomials) the matrix 6^ is usually just the matrix 6^ 
pre-multiplied by a matrix of inner product weights. 

The idea of Balanced POD is to compute the snapshot-based balanced truncation of the system 
in Eq. (11) instead of Eq. (|4]), so that only s adjoint simulations are needed. It is easily shown 
that the systems in Eqs. ( [10] ) and ( [TT] ) have the same observability Gramian using the fact that for 



any projection P, we have = P. Transforming Eq. ( 1 1 ) to balanced coordinates and writing 
X = $ia is obtained as follows: 

a = \I/tA$ia + "^iBu 

(13) 

y = e+C$ia. 

The inverse transformation matrix : M'' ^ A" and the transformation matrix $1 : W ^ X are 
n X r, r being the number of states we want to retain in the system, which we will refer to as the 
rank of the model. Note that r < p, where p is the number of non-zero HSVs. For simplicity, we 
will assume from now on that C = In, i.e. the output of the original system is the full state (this is 
the case in fluid simulations in which we need to know the entire flow field). We can then represent 



the output of Eq. ( 13 ) as y = Qf^iz, which is now the vector of time coefficients of the s standard 
POD modes obtained from the impulse response of the system. For fluids systems the full flow 
field output of the model can be recovered from these coefficients and the corresponding modes. 
For a given dimension of the output projection, all BPOD models will have s outputs regardless of 
the model rank r, while the number of POD model outputs is equal to r at each rank. The effect 
of output projection on model performance will be illustrated in Sec. |lV]and Sec. |Vj 
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III. APPLICATION TO TRANSITIONAL CHANNEL FLOW 



A. Governing equations 



For shear flows, the linearized equations may be conveniently written in terms of the wall- 
normal velocity v and the wall-normal vorticity r] (see, for instance, ll33ll '). The other variables 
(e.g., streamwise and spanwise velocities u and w) may then be computed using the continuity 
equation d^u + dyV + dzW = and the definition of wall-normal vorticity. In these coordinates, 
the linearized (nondimensional) equations have the form 

1 



dt + Ud, - —A 
Re 







-U 



dz 



(14) 
(15) 



Here, Re = UcS/u is the Reynolds number, where u is the kinematic viscosity, 5 is the half- 
width of the channel, and A = + dy + is the Laplacian. Uc is a characteristic velocity, 
which for linearized channel flow is the centerline velocity of the laminar profile U{y). The prime 
indicates differentiation with respect to y. The first equation is the Orr-Sommerfeld equation and 
the second one is known as the Squire equation. It was first shown numerically by Orszag [27] that 
the Orr-Sommerfeld equation for channel flow is stable up to Re ~ 5772, when an exponentially 
unstable eigenmode first arises. The Squire equation has stable eigenmodes for all values of Re. 
Still, complex behavior due to the non-normality exists for stable eigenmodes. The term on the 
right hand side of the Squire equation represents tilting of the spanwise component of the vorticity 
of the mean flow (which here is just U') by the strain rate dv /dz jSl, which gives rise to wall- 
normal vorticity. In the limit of high Reynolds number, the perturbation growth is dominated by 
this process, in particular for streamwise-constant perturbations. While the system also exhibits 
phenomena such as degeneracies and resonances [HI [HI, non-normality has been shown to have a 
dominating effect on the energy growth ll28ll . 

In operator form, we can represent the equations using more compact notation as follows: 



d_ 
dt 



-A 




V 




Los 




V 


/ 




V 




-U'd, LsQ 




V 



(16) 
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Re 



-A' 



where 

Los 

~ - ■ Re- 

are the Orr-Sommerfeld and Squire operators, respectively. If we define 



Ud^A - U"d^ 
1 



A 







-1 




A = 


-A 




Los 








/ 




-U'd, LsQ 



(17) 



with no-slip boundary conditions, we can write the system in standard state-space form: 

X = Ax + Bui + Fu2 
y = Cx 

where B and F represent the spatial distributions of the actuators and disturbances respectively. 



(18) 



with Ui{t) and U2{t) being the corresponding input vectors (the time-dependent amplitudes of the 
columns of B and F). The actuation and the disturbances are equivalent mathematically as they are 
both inputs to the system. We note here that the impulse-state responses are given by Xi (t) = e'^^B 
and X2{t) = e^^F, and the adjoint system impulse-state responses for the full system are given by 
z{t) = e^^*C"^. Therefore, to obtain the POD basis needed for BPOD, we simulate the system 



given by Eq. ( 14) with a given perturbation or actuator as initial condition until the response has 
decayed to negligible levels, so that the matrix XX^ that can be formed from the snapshots will 
closely approximate the controllability Gramian given by Eq. (|5]), where the integral extends to 
infinite time. Of course, computation of the matrix XX^ is intractable for very large systems, so 
we compute POD via the method of snapshots, forming the smaller matrix X^X and following 



the procedure described in Sec. II A 



B. Inner product on the state space 

To determine the corresponding adjoint equations, one first needs to define an inner product on 
the vector space X of flow variables (v, rf). Since balanced truncation is independent of the choice 
of inner product used to define the adjoint (see Appendix |A]), we may choose an inner product 
which is convenient for numerical computations. Let us define the inner product 

(f2,?72))M = / {-viAv2 + r]ir]2)dxdydz, (19) 
Jn 
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where Q denotes the fluid volume. Note that, letting M : X ^ X denote the matrix operator on 



the left hand side of Eq. ( 16 1, this is just the L2 inner product of {vi, rji) with M(t>2, 772)- This 
inner product is different from the standard energy inner product used in analyzing perturbations 
through Fourier decomposition [[3l[8l, as there is no re-scaling at each wavenumber. 

With this definition of the inner product, the adjoint equations are easily found by integration 
by parts: 

(20) 



where 



d 


-A 




V 




Lhs U'd^ 




V 


dt 


/ 




V 




L*sQ 







'OS 



'SQ 



-Ud^A - 2U'd^dy 



Ud, + —A. 
Re 



Re 



-A' 



C. Inner product on the output space 



Although the time evolution of the linearized disturbances is fully described by the wall-normal 
velocity-vorticity formulation, the output of the system can be chosen to be in different variables. 
When using POD, the choice of inner product can have a large impact on the results. If the output 
of our system is only the velocity-vorticity field, the standard L2 inner product can be used. For 
our system, since the other two velocity components can easily be recovered using continuity and 
the definition of vorticity, we can choose the full velocity field to be the output, and use the energy 
inner product given by 



(ui,U2)= / {uiU2 + V1V2 + W1W2) dx dy dz, 
Jq 



(21) 



This choice is more intuitively appealing, since the POD modes for the output projection will 
capture the true kinetic energy of the perturbation. We therefore define the output space y in our 



system as the space M" together with the inner product defined by Eq. (21 ). We note here that the 
space X is also M", though endowed with a different inner product (the M-inner product described 
in the previous section). 
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D. Numerical Methods 



The simulations were performed using a linearized version of a fully nonlinear DNS code 
using the spectral method described by Kim et al. [17|, with periodic boundary conditions in 
the streamwise and span wise directions. The linearized code was verified against the analytic 
time evolution of Orr-Sommerfeld eigenfunctions and optimal perturbations, and the resolution 
of each simulation was checked by varying the time step and grid resolution. The size of the 
computational box was 2n x 2n in the streamwise and spanwise directions for all simulations. 
Standard LAPACK routines were used for the computation of POD and balanced POD modes, 
as well as for the comparison of subspaces. The reduced-order models were integrated using the 
standard fourth-order Runge-Kutta scheme. All computations were done using Fortran 90 and 
MATLAB. A code by S.C.Reddy [|33l was used to compute the initial conditions for the single- 
wavenumber perturbations. The integration weights derived by Hanifi et al BUl were used for the 
computation of inner products on the Chebyshev grid. 



IV. RESULTS: SINGLE- WAVENUMBER PERTURBATIONS 



We start by investigating the system given by Eq. ( 18 1 without actuation and only in the pres- 
ence of disturbances (without the Bui term). In order to validate the numerical methods, we 
first obtain BPOD models from three-dimensional simulations of simple and well-known single- 
wavenumber perturbation cases, described by Butler and Farrell Ol and also investigated by 
Schmid and Henningson [33 1. The general form of such disturbances is given by 

g(x,y,z,t)=g(i/,t)e(^'^^+^'^^), (22) 

where q{y, t) = [v{y, t) r]{y, t)]'^. The standard approach to such perturbations is to compute the 
time evolution of q{y,t), which fully describes the system, since the velocity components u and 
w can easily be computed. For this one-dimensional problem, standard algorithms for computing 
balanced truncation are computationally tractable. Therefore, we are able to compare the models 
resulting from exact balanced truncation (which for the 1-D case can be computed in MATLAB 
using standard algorithms [fT9ll ) to BPOD models obtained from three-dimensional simulations 
of the real part of the full field. Re {q{x,y, z,t)} at a particular wavenumber pair (a,/3) on a 
large grid, similar to the comparison made by Rowley [[30l for a streamwise-constant perturbation. 
We note that for a given wavenumber pair the comparison between BPOD and exact balanced 
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FIG. 1: (a) Kinetic energy growth for the optimal perturbation at wavenumber a = 1, f3 = 1 at Re = 1000. 
(b) The a = 1, P = 1 optimal perturbation, showing streamwise velocity u (complex). 

truncation can be done only using 1-D simulations, but we also performed 3-D simulations in 
order to verify our codes. We also note that, since the outputs of the output-projected system and 



the reduced-order models are coefficients of POD modes, the C matrix in ( [18] ) was modified so 
that the output of the full system is in the POD basis as well. This way, a meaningful comparison 
between the balanced truncation of the full system and BPOD is obtained. 

The initial conditions were computed using the method described by Reddy and Henning- 
son ['281 and their energy growth was verified against values reported in that work. While 
streamwise-constant perturbations exhibit the largest energy growth, three-dimensional perturba- 
tions exhibit more interesting dynamics. (Here, by "three-dimensional," we mean that the per- 
turbations have components in both streamwise and spanwise directions, although the problem 
can still be treated as 1-D in the wall-normal direction as described above.) We focus on the 
a = 1, P = 1 perturbation at Re = 1000, whose energy growth is shown in Fig. [TJ The compu- 
tational grid used in the three-dimensional simulation was 16 x 65 x 16, corresponding to 33280 



states in the system given by Eq. ( 18 1. Balanced truncation of the one-dimensional problem with 
65 Chebyshev modes is easily and accurately computed using the algorithm described in li30ll so 
that BPOD performed on the large system can be compared to exact balanced truncation. 
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5 10 15 5 10 15 

J 3 

(a) (b) 

FIG. 2: (a) The first 15 POD eigenvalues for a = 1,^5 = 1 initial perturbation at Re = 1000. (b) The 
first 15 Hankel singular values (HSVs) for: four-mode (A) and eight-mode (o) output projections and full 
balanced truncation (o) for the same case. 

A. Mode subspaces 

It was found that 500 equally spaced snapshots are sufficient for accurate computation of the 
POD modes, since for a larger number of snapshots with finer spacing there is no considerable 
change in the eigenvalue spectrum or the corresponding modes. We see from Fig. |2] that the 
most significant eigenvalues and the corresponding modes typically come in pairs, representing 
traveling structures that are 90 degrees out of phase. The first pair of modes contains 90.45% of 
the energy, while the first three pairs contain 99.6% of the energy. For the balanced POD models, 
a four-mode and eight-mode output projections were chosen, corresponding to respectively 98.3% 
and 99.9% of total energy contained in the POD modes. 

We also notice that the HSVs (Fig. [2]) come in pairs, indicating that the most significant modes 
in the BPOD mode basis are again traveling structures similar to the standard POD modes. It is 
important to include these pairs of modes in the reduced-order models, as stability of the models 
for balanced truncation is guaranteed only if cr,. > where r is the rank of the model [j4l. While 
for standard POD modes there is no such requirement, mode pairs should always be included in the 
models on physical grounds. We also notice that the number of HSVs for each output projection 
that is equal to the full balanced truncation HSVs is approximately equal to the output projection 
rank. The same observation was made by Rowley [|30l . although there is no proof of this property 
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FIG. 3: Streamwisevelocity for the first POD mode, balancing mode and adjoint mode for the a = 1,P = 1, 
Re = 1000 initial condition. 

at this point. 

The first POD mode is shown in Fig. [3] together with the first balancing and adjoint modes 
from a four-mode output projection. Figure |4] shows the streamwise velocity of the sixth and 
tenth balancing modes, illustrating the effect of the choice of output projection rank. The first four 
balancing modes from BPOD look identical for both output projections, while the sixth mode is not 
very accurately captured by a four-mode output projection. Both output projections do not capture 
very accurately the higher modes such as the tenth mode, which is also illustrated by the HSVs in 
Fig.[2j As we show below, this inaccuracy does not significantly affect model performance, since 
these higher modes are not very significant dynamically. 



In this single-wavenumber case, the exact eigenvalue spectrum of the A matrix from Eq. (17) 
at a given Reynolds number can easily be computed. We note here that the eigenvalues of the 
matrix A and therefore the poles of the corresponding transfer function are independent of the 
initial condition (which is just the B matrix for our impulse response simulations). Figure |5] 
shows the spectra of the full operator and three reduced-order models of different rank for the 
a = 1, P = 1 perturbation. Since the spectra are symmetric about the real axis, we only show 
the upper half of the complex plane. We see that, while the representation of the full spectrum 
improves for both methods as the rank increases, BPOD captures more accurately some of the 
most slowly decaying eigenvalues, which have the most influence on the dynamics of the system. 
For the rank four model, the standard POD model appears to be marginally stable, while the BPOD 
model approximates closely the eigenvalue closest to the origin. At higher order, the standard POD 
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FIG. 4: Streamwise velocity for (a) the sixth balancing mode and (b) the tenth balancing mode for BPOD 
with two different output projections and for full balanced truncation. 
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FIG. 5: Spectrum of the full operator and reduced-order models for rank (a) 4, (b) 8, (c) 30. The BPOD 
modes are from the eight-mode output projection. Symbols: BPOD (□), POD (0)> full operator Only 
the most important part of the full spectrum is shown. 



models improve and capture approximately the same eigenvalues as the BPOD models of the same 
rank. It is also important to notice that some of the eigenvalues of the full system are never captured 
by the reduced-order models. These eigenvalues correspond to uncontrollable eigenmodes of the 
full system, and can never be excited by this particular perturbation. 
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FIG. 6: (a) a = 1, /? = 1 optimal perturbation at Re = 1000, eight-mode output projection, 4-mode and 
8-mode models. Full simulation (+), 4-mode POD (Q), 4-mode BPOD (□), 8-mode POD (o), 8-mode 
BPOD (+) (b) First two outputs, symbols as defined in (a). 

B. Impulse response 



We next compare the impulse response of the system to that of the reduced-order models. The 
impulse response of a linear system is important, since the response of the system to any input can 
be found from the convolution of the impulse response with the input. Figure[6]shows the capturing 
of the growth of kinetic energy by POD and BPOD models, as well as the first two outputs of the 
reduced-order models. The poor performance of POD at low orders for the traveling structure 
perturbation is evident. Even the eight-mode POD model, which captures the energy growth well, 
does not accurately capture the phase of the oscillations. 

Figure |7] shows the 2-norms of the error between the impulse responses of the reduced-order 
models and the full simulation ||G — G^ll, normalized by the 2-norm of the impulse response 
of the full simulation. This figure is a clear demonstration of the effect of output projection. A 
four-mode output projection means that we are effectively performing balanced truncation on the 
dynamics of the first four POD modes of the full system. The dashed lines in the figure indicate the 
2-norms of the error between the full dynamics of the output-projected system and the full system. 
As the rank of the BPOD models is increased, the dashed lines, which are the limit of accuracy, 
are reached fast. As already seen in Fig. [6[ for very low-order models, standard POD is clearly 
outperformed by balanced POD. We see that ten POD modes are needed to match the performance 
of a four-mode output projection BPOD model at the same rank. The slow improvement in POD 
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FIG. 7: Error 2-norms for the a = l,/3 = l,i?e = 1000 perturbation for full balanced truncation, POD and 
BPOD at two output projections. 

model performance indicates that the dynamics of the perturbation can not be represented only 
by retaining the first few most controllable (POD) modes. Adding new BPOD modes beyond 
rank eight and ten (for four- and eight-mode output projections, respectively) does not improve 
the model performance noticeably, since the dynamics of the output-projected system is already 
captured fully. It is also important to note that the performance of the BPOD models is identical 
to that of full balanced truncation almost until the rank at which BPOD model error norms level 
off due to the output projection. This indicates that the higher balancing modes which are not 
computed accurately due to the approximation inherent in the output projection (such as those 
shown in Fig. |4]) do not significantly influence the reduced-order model performance, the main 
hmitation being the capturing of the full system by the output projection. 



C. Frequency Response 

The frequency response encompasses system behavior over the complete range of possible 
forcing, and is perhaps the best indication of overall system performance. Therefore, from the 
control designer's point of view, having a low-order model that represents well the frequency 
response of the original system is of key importance. Frequency response of single-wavenumber 
perturbations was investigated by Schmid and Henningson [|33l using the resolvent norm, where 
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at each frequency the maximum amplification over all initial conditions is computed. Here the 
frequency response of the system with a given actuator or perturbation is of interest. 

A standard way of representing synthesized frequency response for MIMO (multiple input 
multiple output) systems is a plot of the maximum singular value of the transfer function matrix 
max{a{H{iiu))) as a function of frequency, also known as a singular value Bode plot. Fig. [s] 
(a) and (b) shows such plots for the a = 1, P = 1 perturbation and is a clear demonstration 
of the advantages of BPOD for capturing the dynamics of the system. We see that even for a 
two-dimensional model the resonant peak is captured well by the model, while for standard POD 
the peak is very narrow, with very low response at other frequencies. This behavior is typical 
of balanced truncation, as shown in [4J - the first modes to be captured are the ones which are 
most significant dynamically, while the correct response is gradually built up in less significant 
frequency bands as more modes are added. For standard POD models, on the contrary, the re- 
sponse improves incrementally at all frequencies as more modes are added and a higher number of 
modes is needed to accurately capture the resonant peak. For ten mode models, both standard POD 
and BPOD capture well the frequency response (not shown in the figure), with BPOD frequency 
response being almost indistinguishable from the full system one. Standard POD frequency re- 
sponse also includes spurious non-physical peaks at low order, which correspond to eigenvalues 
very close to the imaginary axis for low order of truncation, as seen in Fig. [5] (a). 

We are also interested in the worst-case error between the reduced-order model and the full 
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FIG. 9: Infinity error norms for POD, exact balanced truncation and BPOD with tiie infinity error bounds. 

simulation, which is known as the infinity norm of the system. Balanced truncation has apriori 
error bounds for the infinity norm flU . The error Hoo lower bound for any reduced-order system is 

\\G — Gr\\oo ^ O'r+l, (23) 

where G{s) is the full model, Oris) is the reduced-order model with state dimension r and aj is 
the j-th Hankel singular value (in decreasing order). The upper bound for the error is given by 

IIG-ailoo < 2J:]^,^,aj. (24) 

The upper bound on the error can be very close to the lower bound if the HSVs decrease fast. 
Figure Fig. |9] shows the infinity norm of the error transfer function between the full system and 
the reduced-order model as a function of model rank for the first fifteen orders of truncation. The 



infinity norms for exact balanced truncation lie within the theoretical bounds given by Eqs. (23) 



and (24), while for BPOD, for each of the two output projections, the norms stay within bounds 
up to approximately the rank of the output projection, analogous to the two-norms in Fig. |7| 
The infinity norms for standard POD at low rank are considerably higher than those for balanced 
truncation and BPOD, corresponding to the frequency responses shown in Fig. [81(a) and (b). 
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D. Variation of Reynolds number 



Another very desirable feature of a reduced-order model is good performance for off-design 
values of the system parameters. We would like the models to remain valid for a wide range of the 
model parameters, or at least for the range appropriate for the physical application of the model. 
The only parameter we are considering in our models is the Reynolds number, so the response of 
models was compared to the full simulation when Re is changed. Separating the operators from 



Eq. ( 16 1 into convective and diffusive parts, we can re- write the state-space equation as 



X = AcomX + -^AdifiX + Bu. (25) 
Re 



We can then separately project the matrices A^om and Adiff as in Eq. ( 13 1 at any Reynolds number 
onto the POD and BPOD modes obtained at Re = 1000, the B matrix being just the initial 
condition at Re = 1000. Figure [10] shows the performance of 12-mode standard POD and BPOD 



models when the value of Re in Eq. (25) was changed to 2000 and the impulse response of the 
resulting models was compared to the impulse response of the full system. This rank of the model 
was chosen since both models perform well at the design condition of Re = 1000. We see that the 
BPOD model eigenvalues stay closer to the full simulation eigenvalues (which move as well), and 
also remain in the left half of the complex plane, while for Re = 2000 the standard POD model 
becomes unstable. This indicates a greater range of validity for BPOD models and better stability 
at off-design conditions than standard POD. 

V. RESULTS : THREE-DIMENSIONAL LOCALIZED ACTUATOR 

We next consider a localized body force actuator in the center of the channel which cannot be 



described by a one-dimensional problem. This case corresponds to Eq. (18) without the Fu2 term, 
with the input matrix B representing the velocity field in Fig. [TT] Individual localized disturbances 
to channel flow were investigated by Henningson et al [fTOl . Since balanced truncation involves 
the approximation of the system's Gramians (although in BPOD we do not actually compute the 
Gramians themselves), we are interested in following both the forward and adjoint impulse-state 
responses until all transients have completely decayed. The computational box necessary for fol- 
lowing individual localized disturbances long enough in time would be very large, and we instead 
consider a periodic array of small disturbances in the channel center. It should be noted that the 
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FIG. 10: Left: The comparison of spectra of the full operator at a = 1, /? = 1 to the spectra of rank 12 
models as the value of Re is increased to 2000. Right: The performance of corresponding reduced-order 
models at Re = 2000 - first two outputs. See text for detailed description. Symbols: full (x), BPOD (□), 
POD (O). 

behavior of this periodic array can be quite different from the behavior of a single localized distur- 
bance, in particular considering the energy growth, since the periodic array quickly develops into 
a streamwise-constant vortex. The exact form of the initial condition considered here is 

^2 ■ 



v{x,y,z,0) = A 1 - — e^-' -y /<\cos{ny) + 1) 



(26) 



' + 



where Xc,0,Zc are the coordinates of the center of the computational domain and 
(z — ZcY. The wall-normal vorticity is zero. This form was picked in order to satisfy the condition 
that the mean perturbation velocity is zero at each wavenumber. The (cos(7ry) -|- 1) term was added 
to make the field satisfy exactly the boundary conditions v{±l) = Vy(±l) = 0. The amplitude A 
was set to 1 for this simulation, and the parameters a and ay were set to a = 0.7 and ay = 0.6. 
The Reynolds number chosen for this simulation was Re = 2000. The traveling structure rapidly 
develops into a streamwise-constant structure, since the growth of wall-normal vorticity results in 



the development of streamwise streaks (see Fig.[TT]). 

The grid size was 32 x 65 x 32, corresponding to 133,120 states for the full {v, rj) system. The 
simulation was ran for 1200 dimensionless time units (t = f^Uc/S), and the timestep used was 
At = 0.004. During this time, the energy of the initial disturbance decayed to approximately 1.5 
percent of its initial value. The POD modes were taken over 1000 snapshots, with fine spacing 
between snapshots for the initial period in order to capture the traveling structures well and larger 
spacing once the streamwise structure was developed, after it was verified that POD eigenvalues 
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FIG. 1 1 : The development of the wall-normal velocity of the perturbation given by Eq. ( 26 1 at t = (left), 
t = 14 (middle) and t = 160 (right) which corresponds to the maximum energy growth. Positive velocity 
is light, and negative velocity is dark. 
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(a) (b) 
FIG. 12: (a) The first 20 POD eigenvalues for the Gaussian- like disturbance impulse response, (b) The first 
20 HSVs for five-mode (o) and ten-mode (A) output projections. 



and the corresponding modes do not change significantly if more snapshots are used. Fig. 12 (a) 
shows the POD eigenvalues of the impulse response. The first five modes contain 99.72% of the 
perturbation energy, and the first ten modes contain 99.9% of the energy. In this case the spectrum 
contains both streamwise-constant (and nearly- streamwise constant modes) as well as traveling 
structures due to the initial transient. The first three modes are streamwise-constant structures, 
while the fourth and the fifth modes correspond to a traveling structure, which accounts for only 



0.40% of the total energy. Modes one, four and five are shown in Fig. 13 

Next, the adjoint simulations were ran and the BPOD procedure was performed on a five-mode 
output projection, containing only the most important traveling structure, as well as on a ten-mode 
output projection. These ranks were chosen due to large drops in energy significance after the 



fifth and tenth mode, as shown in Fig. [12^ (a). Fig. 12 (b) shows the HSVs for these two output 
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FIG. 13: The first, fourth and fifth POD modes for the locahzed actuator, showing streamwise velocity. The 
iso-surfaces show half of the maximum (light) and minimum (dark) value. 




FIG. 14: Top row: primal modes one, four and five from balanced POD, showing streamwise velocity for 
the localized perturbation. The modes are from a five-mode output projection. The iso-surfaces show half 
of the maximum (light) and minimum (dark) value. Bottom row: the corresponding adjoint modes. Note 



the similarity between the primal modes and the corresponding POD modes in Fig. 13 



projections. We notice that the HSVs are equal for the pairs of modes 4-5 and 7-8 for five-mode 
output projection, corresponding to traveling structures in the basis of BPOD modes. Even more 
interestingly, for the ten-mode output projection, HSVs for the modes 4-6 are equal. Although the 
stability of balanced truncation models is guaranteed only when ar+i < (Xr, where r is the number 
of states retained [4], 4-mode, 5-mode and 7-mode models for both output projections were found 
to be stable. The model error for impulse response, however, decreases significantly if both modes 



corresponding to a traveling structure are included, as will be shown in Sec. VA (see Fig. 16) 



Balanced POD modes one, four, and five, and the corresponding adjoint modes for the five-mode 

is 



output projection are shown in Fig. [141 Note that the structure of modes four and five in Fig. 14 



almost identical, except for a spatial phase shift of exactly one quarter of the periodic domain. 
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A. Impulse Response 



Figure[T5](a) shows the perturbation energy growth as captured by three different standard POD 



models. It was observed that the inclusion of modes which come in pairs (see Fig. 12 (a)) in the 
basis used to form the reduced order models does not change the system behavior appreciably - 
the response of a model including modes 1-9 (not shown in figure) is virtually indistinguishable 
from the response of the model including only the first three modes. Hence, the traveling structure 
modes do not contribute significantly to the dynamics of this perturbation. The inclusion of the 
tenth mode, which is streamwise-constant, improves the performance significantly, and the model 
composed of only the first three modes and the tenth mode performs as well as one including 
the first ten modes. In the same fashion, including the mode pairs 11-12, 13-14 and 15-16 does 
not affect the model performance. Including the seventeenth mode, which is also a streamwise- 
constant mode, improves the performance further. The tenth and the seventeenth mode correspond 
to 0.025% and 0.0074% of the total energy. The low-order standard POD models were found to 
capture poorly the initial condition of the full simulation (this will be discussed in more detail 



in Sec. VD), so they were also started from different initial conditions at later times (before or 
around the peak energy growth), when the projection of the simulation onto POD modes is close 
to the full simulation data, and they still did not capture the correct peak and the subsequent decay 
of the energy. 

On the other hand, the performance of very low-order BPOD models is significantly better. Fig- 
ure [15] (b) shows the perturbation energy growth as captured by three different models. Although 



the two-mode model does not accurately capture the initial condition, it does represent the energy 
growth at later stages reasonably well. A three-mode BPOD model captures the kinetic energy 
of the full simulation very well except for the initial period. While more modes are needed to 
capture exactly the initial transient, if only the energy growth is of interest, the three-mode model 
is sufficient. This striking difference is an illustration of the advantage of balanced truncation - for 
standard POD it is difficult to know apriori which modes will be important for the system dynam- 
ics as demonstrated above and a good low-order model was found only after a careful examination 
of the mode basis which provided some insight into the underlying physics. 



Figure 16 shows the error HG — G^H for BPOD using the two output projections and for stan- 
dard POD. The error for the BPOD models quickly reaches the asymptotic values dictated by the 
output projection, although more modes are needed compared to the optimal perturbation case de- 
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FIG. 15: (a) Three, four and five-mode POD models formed from the indicated modes (b) Two-mode and 
three-mode BPOD models. The very low-order BPOD models do not capture very well the initial transient, 
as shown in the inset. The BPOD models are from the ten-mode output projection. 




FIG. 16: Error 2-norms for localized actuator, showing POD models and BPOD at two output projections. 

scribed in the previous section due to the more complex dynamics. Standard POD starts to match 
the performance of the ten-mode output projection BPOD only around rank 30, and also varies a 
lot with the model rank. This corresponds to the already observed fact that dynamically impor- 
tant POD modes are not highly ranked in terms of energy. Whenever the POD modes come in 
pairs, including only one of the modes results in deterioration of model performance. Eventually 
standard POD has better performance than BPOD, however recall that these POD models have r 
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outputs while the BPOD models have only s outputs. It should also be noted that some POD mod- 
els exhibit sustained or very slowly decaying oscillations, and that the corresponding two-norms 



are in fact infinite. Since the simulation time for Fig. 16 is finite, the two-norms of such models 
appear to be large but finite as well. Although the four-mode BPOD models are stable for both 
output projections, including just one of the modes corresponding to a pair of equal HSVs 4-5 
for the five-mode output projection deteriorates model performance, while we see a large decrease 
in the error when the fifth mode is included, as well as when we include subsequent pairs. For 
the ten-mode output projection, there are three equal HSVs 4-6, and a significant decrease in the 
error is seen only when we include all three of those modes (in particular, the error norm of the 
five-mode BPOD model is significantly larger than that of the four- mode model). 



B. Frequency Response 



Figure 17 shows the singular value Bode plots of standard POD and ten-mode output projec- 
tion BPOD models for the localized disturbance. The frequency response of the 50-mode BPOD 
model, which is a very close approximation of the frequency response of the actual disturbance, 
has the shape of a low-pass filter with a break frequency of 0.01 rad/s with two resonant peaks 
near 1 rad/s, which are similar to the peak observed for the single- wavenumber traveling struc- 
ture perturbation in the previous section. We see that standard POD models again have spurious 
peaks at low model ranks. The addition of mode pairs corresponding to the traveling structures is 
necessary in order to reproduce the peaks around 1 rad/s for both standard POD and BPOD, how- 
ever BPOD captures those peaks with only the triple 4-6 and the mode pair 7-8, as well as modes 



9 and 10 (Figure 17 (b)) while all standard POD modes 1-17 are needed to reproduce the same 
peaks and there are still spurious peaks. Since the peaks correspond to the low-energy traveling 
structures, it is not surprising that only a three-mode BPOD model performs so well in capturing 
the kinetic energy of the full simulation, as shown in the previous section. On the other hand, 
if the frequency response of the actuator around the frequency of 1 rad/s needs to be captured 
accurately, the higher BPOD modes need to be included. 
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FIG. 17: Singular value Bode plots for POD (a) and BPOD (b) models for the localized disturbance. The 
models are compared to a 50-mode BPOD model. The close-up in (b) shows that a six-mode BPOD model 
is needed to capture the larger resonant peak, and a ten-mode BPOD model captures both peaks. 

C. Variation of Reynolds number 



Figure 18 shows some of the eigenvalues of the 17-mode POD and BPOD models as the 



Reynolds number is increased. As in Section IV D we use the modes from the design condi 



tion (Re = 2000 in this case) and form the models using Eq. (25 1. Both POD and BPOD models 
have eigenvalues on the real axis very close to the origin, which remain stable and correspond 
to the slow evolution of the streamwise-constant structures. At each Re, the eigenvalues of both 
models move towards the right half of the complex plane and while the BPOD model always 
remains stable, the standard POD model first appears marginally stable at Re = 2500 and then 
unstable at Re = 3000. The effect of the eigenvalues that move to the right half of the complex 



plane is clearly seen in Fig. 18 (b). A model that includes modes 1-17 grows unstable quickly 
at Re = 3000, showing that inclusion of modes which at design condition do not contribute sig- 
nifcantly to the overall dynamics can significantly deteriorate the performance of the model at 
off-design condition. This can also be seen from the frequency responses shown in the previous 
section - even at design condition, the spurious high peaks correspond to marginally stable modes. 
Although stable, the 1,2,3,10,17 POD model is highly inaccurate at Re = 3000, with large peaks 
in the kinetic energy which decay very slowly, indicating the high sensitivity of those POD mod- 
els which remain stable to a change in the Reynolds number. On the other hand, the three-mode 
BPOD model is still remarkably close to the full system. 
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FIG. 18: (a) The eigenvalues of 17-mode models at Re = 2000, Re = 2500, and Re = 3000. (b) The 
performance of the models at the off-design condition of Re = 3000. 

It is important to note here that as the Reynolds number is increased, the nonlinearity will have 
a stronger effect on the development of the disturbance and the reduced-order model may not be 
valid for higher Re in the first place. The comparison of the linear perturbation growth with a full 
nonlinear DNS solution is essential for a true validation of the models for controls applications, 
since we may be modeling the linearized flow well, but the linearized flow may not be a good 
approximation to the actual flow. This comparison is subject of current work. 



D. Capturing of actuation 

An important property of a reduced-order model is how well it captures the effects of the actu- 
ator in the original system, especially for models that are intended for developing controllers. In 
order for a reduced-order model to capture the effect of an actuator, it is necessary at a minimum 
for the input term in the equations (Bu in Eq. (|4])) to be contained in the subspace used for pro- 
jecting the equations. Note that here, even for the standard POD case, the effect of the "actuator" 
is partially included, since the dataset used for POD is generated by an impulsive input. One way 
to measure the degree to which the input "directions" are captured by the modes used in the model 
is to compute the projection of the columns of the input matrix B in Eq. (|4]) onto the basis modes. 
In the system we are considering here, i? is a single column vector, representing the initial distur- 



bance given to the system (or actuation via a body force in the center of the domain). Fig. 19 shows 
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the norm of the projection ||P,.i?||/||i?|| of the standard and BPOD modes onto the input vector 
B, which is just the initial condition for each simulation. The balancing modes clearly capture 
the input direction with many fewer modes than standard POD: even very low-order models have 
a significant norm after projection, and in fact the norm of B after projection is almost always 



greater than the norm of B due to the non-orthogonal projection, as shown in Fig. 19 Any orthog- 
onal projection such as P must satisfy ||-Px|| < while for a non-orthogonal projection we 
may have > which is the case for the first several BPOD modes. Clearly, the B matrix 

has a very small projection onto the POD modes for standard POD unless many modes are taken, 
so it is impossible for very low-order POD models to capture the response of an actuator without 
introducing more modes (such as the B matrix itself, Krylov subspaces, or shift modes [jlSl ). 




(a) (b) 
FIG. 19: Norm of the projection of the B matrix (a single column vector) onto subspaces used for reduced- 
order models, (a) a = 1, /? = 1, (b) localized disturbance. The diagram in (a) illustrates the non-orthogonal 
projection used in BPOD. 



E. Subspace comparison 

The BPOD procedure uses both a different projection and a different set of modes in order 
to form reduced-order models, and we next look at a comparison of the two mode subspaces. 
A way to compare two subspaces is to compute Tr^PAPsPA) = T, where Pa and Pb are the 
corresponding projection operators [[51. The trace of the matrix Tr^PAPsPA) as a function of 



the subspace rank is shown in Fig. 20 for the five-mode and ten-mode output projections where 



Pa and Pb are the orthogonal projectors onto the POD and BPOD subspaces respectively. The 
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FIG. 20: Plots of Tr{PAPBPA) as a function of the subspace rank r for the locahzed perturbation for both 
output projections. 

value of the trace T is the same as the subspace rank r at low order, indicating very similar modes 



structures, and including both modes from the pair brings the value of s to almost exactly 5. It is 
interesting to observe that for the five-mode output projection, r = T exactly at r = 5, while for 
the ten-mode output projection r = T at r = 10, and that above those values the value of the trace 
is lower than the rank. This can be explained by the fact that BPOD is attempting to approximate 
the output projection of the data of the given rank. It is interesting to note that the subspaces 
including the first three POD and BPOD modes are virtually identical, indicating that the non- 
orthogonal Petrov-Galerkin projection via adjoint modes makes the enormous difference that we 
have seen in the performance of the corresponding models. As mentioned above, the POD basis 
is the basis of the most controllable modes, and is indeed optimal in capturing a given dataset, but 
as we have shown, it can fail to capture the dynamics correctly. 

VI. CONCLUSIONS AND FURTHER WORK 

We have shown that the balanced POD procedure offers several advantages over the standard 
POD technique for modeling linearized channel flow. Balanced POD makes the benefits of the 
systematic model reduction via balanced truncation available for large systems, resulting in mod- 



(see Figs. 13 and [14]). For both POD and five-mode OP BPOD, modes four and five are a pair of 
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els that capture the impulse response, frequency response and the effects of an actuator much 
better than their corresponding POD models. The differences between POD and BPOD are espe- 
cially striking for localized actuators or disturbances which are often of interest in applications. 
As BPOD allows model reduction for very large systems, the localized perturbation models were 
obtained using three-dimensional fields in physical space, without the standard wavenumber de- 
coupling which can be unwieldy for localized disturbances, since model reduction would have 
to be performed at each wavenumber pair resulting in models which are still of relatively high 
rank, and each wavenumber pair would need to be controlled separately. More importantly, this 
approach allows the extraction of dynamically dominant structures in physical space which may 
have contributions from many wavenumbers. This feature of the method will be very important in 
further applications of this method, and in particular in spatially developing flows. 

There is an indication that the non-orthogonal projection onto the balancing modes with the use 
of adjoint modes plays a key role in the model performance. The subspaces spanned by the BPOD 
modes are very close to those spanned by the POD modes (which are optimal in representing the 
impulse response dataset). However, the dynamical models produced by BPOD are quite different, 
because of the non- orthogonal projection (or, equivalently, an orthogonal projection with a new 
inner product, defined by the observability Gramian [|30ll ). 

This work has focused on the linearized system in order to characterize the BPOD model re- 
duction method. A true test of any controller is the application to the full nonlinear simulation, 
and verification of reduced perturbation growth, and this is a subject of our ongoing work. In addi- 
tion, several techniques are available for obtaining nonlinear reduced-order models. One method 
for snapshot-based balancing for nonlinear systems has been introduced by Lall et al [fTSl, but 
this involves considerably more computation than the present method, and is not feasible for large 
systems. A simpler, but presumably less effective approach is to project the full Navier-Stokes 
equations onto the BPOD modes computed for a linearized system; such a procedure involves no 
additional computational expense over the methods presented in this paper. 

Although the body force actuation we discuss here would be desirable in applications, it is not 
as practically feasible as wall blowing and suction, which has been studied extensively and has 
shown promise both in computations and experiment lfT2ll2Tll . We also note that blowing and suc- 
tion could be incorporated into the control framework described here in a relatively straightforward 
fashion, via the lifting approach described in Ref. [[T2l . Finally, the BPOD procedure is applicable 
to other flows that can be linearized, such as Couette flow, boundary layer flows, or flows over 
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airfoils. Moreover, the method can be appUed to any large linear system and can therefore be a 
useful addition to the tools of modem control theory. 
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APPENDIX A: SNAPSHOT-BASED BALANCED TRUNCATION USING A CONTINUOUS AD- 
JOINT 

When computing the exact balanced truncation, the balancing transformation is found from the 
eigenvalue problem W^WqT = TTi^ where Wc and Wo are the controllability and observability 
Gramians with — . We show here that, although the product WcWo does not depend on the 
inner product on the state space used to define the adjoint system, the appropriate weight M needs 
to be included in the computation via the method of snapshots. 

We can represent the weighted inner product of two vectors qi and q2 as 



where the domain of integration Q, is the Hilbert space itself. The star denotes the complex conju- 
gate transpose. The inner product weight M is part of the definition of the Hilbert space itself. We 
define the so-called continuous adjoint of an operator A with respect to this inner product as 



We use the symbol + in order to distinguish the adjoint from the standard matrix transpose A^ . 
From this definition it is easily shown that 




(Al) 



M 



(A2) 



{Ax, z)^ = (x, A^z)^ ^ ^+ - M-^A^M 
{Bu,x)j^^ {u,B+x) B^^B'^M 



(A3) 



(Cx, y) = {x, C+y)^ ^ C+ = M-'C^ 
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In the above, we have assumed that the input and output spaces use the standard (unweighted) 
inner product. Next, we obtain for the Gramians: 



oo POO 



Jo Jo /A4^ 

poo POO V 

Go= e^^'C+Ce'^'dt= / M-^e^^'MM-^C^'Ce^'dt 
Jo Jo 

where Gc and Go denote the Gramians obtained with the weighted inner product. Since the matri- 
ces M and are constant, we can take them out of the integrals, obtaining 

G,Go = WMo (A5) 

Thus, we have shown that balanced truncation does not depend on the choice of the inner 
product used to derive the adjoint system, and this allows us to use a convenient inner product. 
(In numerical simulations the 'simple' discrete adjoint = may in fact be more difficult to 
compute than a continuous adjoint which may retain a similar form of the equations; for instance, 
this is the case for linearized channel flow). 

Next, we consider the computation of balancing and adjoint modes via the method of snapshots. 
From the definition of the empirical Gramians (Eq. ([v])) it is easily shown that = Y'^M 
(recall that the snapshots of the adjoint simulations, which are the columns of Y, are given by 
z{t) = e^^*C+). Thus, we can write the SVD in Eq. § as 



Y^MX = UTy^ (A6) 

If we define the inverse of the balancing transformation as = T.^^^'^U'^Y^ we can easily 
compute the adjoint modes just from the SVD and from the adjoint snapshots. Recall that the 
columns of ^ give the adjoint modes. The two sets of modes will now be bi-orthogonal with 
respect to the M inner product, so that ^^A'/$ = /. 

An alternative, more intuitive explanation is that, since both the direct and the adjoint snapshots 
'live' in the state space, the correct inner product is that including the weight M (which is a part of 
the definition of the Hilbert space in which they reside). It is therefore this weighted inner product 
that should be used for forming the matrix for the SVD. Furthermore, the balancing and adjoint 
modes are bi-orthogonal with respect to this weighted inner product. 
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